// Get stopping positions of muons inside the target
#include "common.h"
#include "../rootlibs/TTrackerHit.h"
#include "../rootlibs/TEvent.h"

void GetOrgMomentum(int runNo, TH1F *h);
void ProcessRun(int runNo, TH2D *h)
{
	printf("Run %d:\n",runNo);
	char filename[20];
	sprintf(filename,"run%05d.root",runNo);
	TFile *file = TFile::Open(filename);
	TTree *t = (TTree*)file->Get("t");
	TEvent* evt = new TEvent();
	t->SetBranchAddress("evt",&evt);
	int nevt = t->GetEntries();
	// Real work follows here
	h->SetXTitle("Muon momentum (MeV/c)");
	h->SetXTitle("X (mm)");
	h->SetYTitle("Y (mm)");
	timeval t1,t2;
	int nstop = 0;
	int nstop0 = 0;
	int nhit_in_target = 0;
	int nstop_in_target = 0;

	gettimeofday(&t1,NULL);
	for (int i = 0; i < nevt; i++) 
	{
		t->GetEntry(i);
		int nhits = evt->GetNoOfHits();
		for (int j = 0; j < nhits; j++) 
		{
			TTrackerHit aHit = evt->GetHit(j);
			bool is_mu = strcmp(evt->GetHit(j).GetPName(),"mu-") == 0;
			if (is_mu && aHit.GetDetID()==0 && aHit.GetHitType()==0)
				nstop0++;
			if (is_mu && aHit.GetDetID()==1)
			{
				if (aHit.GetHitType()==0)
					nstop++;
				if (aHit.GetHitType()==1)
				{
					double pos[3];
					aHit.GetPos(pos);
					h->Fill(pos[0]/1000,pos[1]/1000);
				}
				double pos[3];
				aHit.GetPos(pos);
				bool intarget = pos[0]>-25000 && pos[0]<25000 && pos[1]>-25000 
					&& pos[1]<25000;
				if (intarget && aHit.GetHitType()==1)
					nhit_in_target++;
				if (intarget && aHit.GetHitType()==0)
					nstop_in_target++;
			}	
		}
	}
	gettimeofday(&t2,NULL);
	double timelapse;
	timelapse = (t2.tv_sec - t1.tv_sec)*1000000 + (t2.tv_usec - t1.tv_usec);

	std::cout<<"Stop Plastic: "<<double(nstop0)/nevt*100<<std::endl;
	std::cout<<"Stop Al: "<<double(nstop)/nevt*100<<std::endl;
	std::cout<<"hit inside target: "<<double(nhit_in_target)/nevt*100<<std::endl;
	std::cout<<"Stop inside target: "<<double(nstop_in_target)/nevt*100<<std::endl;
	std::cout<<"Time lapse: "<<timelapse<<std::endl;
}

int main (int argc, char *argv[])
{
	TApplication theApp("App",&argc,argv); 
	SetStyle();

	TH2D *h_mu_hit_target[9];
	for (int i = 0; i < 9; i++) 
	{
		char h_str[20];
		sprintf(h_str, "h%02d",i);
		char h_title[100];
		int p0 = 30 - i*1;
		sprintf(h_title,"Muon hit on target, 0.5 mm scint, 0.2 mm mylar, init: %d MeV/c",p0);
		h_mu_hit_target[i] = new TH2D(h_str,h_title,120,-50,50,120,-50,50);;
		h_mu_hit_target[i]->SetLineColor(1+i);
	}
	
	//GetOrgMomentum(1,h_mu_org);

	for (int i = 0; i < 9; i++) 
	{
		//ProcessRun(1+i, h_mu_hit_target[i]); // 0.5 mm
		ProcessRun(23+i, h_mu_hit_target[i]); // 1.0 mm
	}

	char canvas_str[50];
	sprintf(canvas_str,"Muon hit pos on target");

	TCanvas *c2 = new TCanvas("c2",canvas_str,900,1000);
	//TLegend *leg = new TLegend(0.15,0.6,0.31,0.89);
	//leg->SetFillStyle(0);
	//leg->Draw();
	c2->Divide(3,3);
	for (int i = 0; i < 9; i++) 
	{
		c2->cd(i+1);
		h_mu_hit_target[i]->Draw("colz");
		TLine *l = new TLine();
		l->SetLineColor(kRed);
		l->DrawLine(-25,-25,-25,25);
		l->DrawLine(-25,25,25,25);
		l->DrawLine(+25,25,25,-25);
		l->DrawLine(25,-25,-25,-25);
	}

	char pdffile[50];
	sprintf(pdffile,"pdf/mu_hit_pos_200um_target.eps");
	c2->Print(pdffile);
	//sprintf(pdffile,"pdf/Mu_momentum_indi.eps");
	//c2->Print(pdffile);

	std::cout<<"Done\n";


	// To keep the ROOT windows open after execution
	theApp.Run(); 
	return 0;
}

void GetOrgMomentum(int runNo, TH1F *h)
{
	printf("Run %d:\n",runNo);
	char filename[20];
	sprintf(filename,"run%05d.root",runNo);
	TFile *file = TFile::Open(filename);
	TTree *t = (TTree*)file->Get("t");
	TEvent* evt = new TEvent();
	t->SetBranchAddress("evt",&evt);
	int nevt = t->GetEntries();
	h->SetXTitle("Muon momentum (MeV/c)");

	for (int i = 0; i < nevt; i++) 
	{
		t->GetEntry(i);
		int nhits = evt->GetNoOfHits();
		for (int j = 0; j < nhits; j++) 
		{
			TTrackerHit aHit = evt->GetHit(j);
			bool is_mu = strcmp(evt->GetHit(j).GetPName(),"mu-") == 0;
			if (is_mu && aHit.GetDetID()==0 && aHit.GetHitType()==1)
			{
				double ke = aHit.GetKE();
				h->Fill(GetMuMomemtum(ke));
			}
		}
	}
}
